//
//      +-----^-----+
//      |           |
//      |           | 
//      |     *     |
//      |   (i,j+1) |
//      |           |
//      +-----^-----+
//      |   (i,j)   |
//      |           | 
//      |     *     |
//      |   (i,j)   |
//      |           |
//      +-----^-----+
//                       
//                  |           |           |           |
//                -->-----o----->-----o----->-----o----->-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |  (i,j+1)  |     :     |
//                  |     ^     |    v_N    |     ^     |   
//                  |     :     |     :     |     :     |
//                  |     :     |     :     |     :     |
//                -->-----o---- 3 -- v_n -- 4 ----o----->--  4 = u(i  , j+1)
//                  |     :     |     :     |     :     |    3 = u(i-1, j+1)
//                  |     :     |     :     |     :     |    2 = u(i  , j  )
//                  |    v_W   u_w   v_P   u_e   v_E    |    1 = u(i-1, j  )
//                  |  (i-1,j)  |   (i,j)   |  (i+1,j)  |
//                  |     :     |     :     |     :     |
//                -->-----o---- 1 -- v_s -- 2 ----o----->-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |     :     |     :     |
//                  |     ^     |    v_S    |     ^     |   
//                  |     :     |  (i,j-1)  |     :     |
//                  |     :     |     :     |     :     |
//                -->-----o----->-----o----->-----o----->--
//                  |           |           |           | 
//                  
//               1           3                       2         4 
//   u_w = ( u(i-1,j) + u(i-1,j+1) ) / 2   u_e = ( u(i,j) + u(i,j+1) ) / 2
//   v_n = ( v(i,j) + v(i,j+1) ) / 2     v_s = ( v(i,j) + v(i,j-1) ) / 2
//               
//---------------------  3D  ---------------------
//

namespace Tuna {

template<class T_number, int Dim>
inline
bool CDS_YLES<T_number, Dim>::calcCoefficients(const ScalarField &nut) { 
    T_number dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    T_number dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    T_number ce, cw, cn, cs, cf, cb;
    T_number nutinter;
    T_number dxyz_dt = dx * dy * dz / dt;
    T_number RaGaVol = Rayleigh * Gamma * 0.5 * dx * dy * dz;
   
    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = ( u(i,j,k) + u(i,j+1,k) ) * 0.5 * dyz;
		cw = ( u(i-1,j,k) + u(i-1,j+1,k) ) * 0.5 * dyz;
		cn = ( v(i,j,k) + v(i,j+1,k) ) * 0.5 * dxz;
		cs = ( v(i,j,k) + v(i,j-1,k) ) * 0.5 * dxz;
		cf = ( w(i,j,k) + w(i,j,k+1) ) * 0.5 * dxy;
		cb = ( w(i-1,j,k) + w(i-1,j,k+1) ) * 0.5 * dxy;

//
// nut is calculated on center of volumes, therefore, nut
// must be staggered in y direction:	    
		nutinter = 0.5 * ( nut(i,j,k) + nut(i,j+1,k) );

		aE (i,j,k) = (Gamma + nutinter) * dyz_dx - ce * 0.5;
		aW (i,j,k) = (Gamma + nutinter) * dyz_dx + cw * 0.5;
		aN (i,j,k) = 2 * (Gamma + nutinter) * dxz_dy - cn * 0.5;
		aS (i,j,k) = 2 * (Gamma + nutinter) * dxz_dy + cs * 0.5;
		aF (i,j,k) = (Gamma + nutinter) * dxy_dz - cf * 0.5;
		aB (i,j,k) = (Gamma + nutinter) * dxy_dz + cb * 0.5;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) +
		             aN (i,j,k) + aS (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) +
		             dxyz_dt;	    
//		aP (i,j,k) /= alpha;  // under-relaxation
//		+ (ce - cw)  + (cn - cs) + (cf - cb);	    
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.

		sp (i,j,k) = v(i,j,k) * dxyz_dt - 
		    ( p(i,j+1,k) - p(i,j,k) ) * dxz +
		    RaGaVol * ( T(i,j,k) + T(i,j+1,k) ) +
		    nutinter * ( (u(i,j+1,k) - u(i,j,k) - 
				  u(i-1,j+1,k) + u(i-1,j,k)) * dz +
				 (w(i,j+1,k) - w(i,j,k) - 
				  w(i,j+1,k-1) + w(i,j,k-1)) * dx );
		
//		    v(i,j,k) * (1-alpha) * aP(i,j,k)/alpha; // under-relaxation
	}    
    calc_dv_3D();
    applyBoundaryConditions3D();

    return 1;
}

} // Tuna namespace















